Gridded data in ‘sp’

Again, we can create objects,here start from grid topology definition

topology <- GridTopology(cellcentre.offset = c(1,1,2), cellsize=c(1,1,1), cells.dim = c(3,4,6)) #cellcentre.offset  numeric; vector smin coordinates of cell centre in each dim 
sp_grid <- SpatialGrid(topology)
class(sp_grid)

[1] “SpatialGrid” attr(,“package”) [1] “sp”

str(sp_grid)

Formal class ‘SpatialGrid’ [package “sp”] with 3 slots ..@ grid :Formal class ‘GridTopology’ [package “sp”] with 3 slots .. .. ..@ cellcentre.offset: num [1:3] 1 1 2 .. .. ..@ cellsize : num [1:3] 1 1 1 .. .. ..@ cells.dim : int [1:3] 3 4 6 ..@ bbox : num [1:3, 1:2] 0,5 0,5 1,5 3,5 4,5 7,5 .. ..- attr(*, “dimnames”)=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] “min” “max” ..@ proj4string:Formal class ‘CRS’ [package “sp”] with 1 slot .. .. ..@ projargs: chr NA

What we have made so far

Create Grid object from points

Start to build objects from grid topology, coordinates or read gridded data

points_grid <-  expand.grid(x = 1:3, y = 1:4) 
points_grid 

   x y
1  1 1
2  2 1
3  3 1
4  1 2
5  2 2
6  3 2
7  1 3
8  2 3
9  3 3
10 1 4
11 2 4
12 3 4

sp_points_grid <- SpatialPoints(points_grid)
class(sp_points_grid)

[1] "SpatialPoints"
attr(,"package")
[1] "sp"

Visualise prepared grid

Continuation to full grid

#from SpatialPoints SpatialPixels
pix_grid <- SpatialPixels(sp_points_grid)
#or this way object gridded (SpatialPoints into SpatialPixels)
class(sp_points_grid)

[1] "SpatialPoints"
attr(,"package")
[1] "sp"

gridded(sp_points_grid) <- T
#is objet full grid?
fullgrid(pix_grid)  

[1] FALSE

Continuation to full grid

fullgrid(pix_grid) <- T #set for TRUE and create SpatialGrid
#add atributes (data.frame), to spatial components
atr_grid <- data.frame(matrix(1:12)) #data.frame
atr_grid[1:3,]#first 3 rows of data.frame

[1] 1 2 3

class(atr_grid)

[1] "data.frame"

my_grid <- SpatialGridDataFrame(pix_grid,  atr_grid)
class(my_grid)

[1] "SpatialGridDataFrame"
attr(,"package")
[1] "sp"

visualise prepared grid

Visaualisation with lattice graphics

Import gridded data with sp package

Common way of importing gridded data; calculate summary statistics…

dem <- readGDAL ("dem_1k.asc")###!!!20k

dem_1k.asc has GDAL driver FITS and has 21 rows and 21 columns

#summary statistics
summary(dem@data$band1)

Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 

0,000 6,419 157,996 301,073 511,735 1255,097

With all types of spatial data, statistical procedures from R

As example plot histogram of elevation values

Visualise gridded data with attributes

Export as …

ESRI asc, SAGA sgrd…

writeGDAL(dem, fname="dem.asc", drivername="AAIGrid",  mvFlag=99999)
writeGDAL (dem, fname= "dem.sdat",  drivername="SAGA") 
writeGDAL (dem, fname="dem", drivername="GSBG")

Gridded data with raster package

• Extremely popular for analysis of big rasters

• Data does not need to be in memory

• Easy shift to ‘sp’ classes

• Extract wanted spatial extent or aggregate by polygon, select on points…

Get data from Worldclim database - raster object

library(raster)
x <- getData('worldclim', var='tmin', res=10, lon=5, lat=45)

Plot raster object

library(raster)
plot(getData('worldclim', var='tmin', res=10, lon=5, lat=45))

Extract values for polygons

Values of Digital elevation model for each county

dem_1k.asc has GDAL driver FITS and has 21 rows and 21 columns

dem_r <- raster(dem)
#extracted values for each 1 km cell in polygons
#object of class list, length 21, as many counties
dem_county <- extract(dem_r, HRV_adm1_ETRS)
dem_county[1:20]

[[1]] [1] 722,6663 261,4872 370,5971 773,7128 121,7043 307,7529 601,9179 259,0665

[[2]] [1] 157,0298 122,3399 136,0429 153,9509 157,9960 149,9645 356,1018

[[3]] [1] 125,55239 143,63921 144,06564 97,78848

[[4]] [1] 98,8438 387,1294

[[5]] [1] 196,9077 168,1451 122,0084

[[6]] [1] 185,08521 377,42798 163,93018 281,47314 64,61814 173,73517 24,62883

[[7]] [1] 145,2687 143,5296 336,7780 201,0698 206,9406 757,8558 546,8431 [8] 352,1976 266,6391 691,3894

[[8]] [1] 198,0352 150,9576 123,2260 138,8395 153,0450

[[9]] [1] 285,3583 284,5465 202,5881 237,4919

[[10]] [1] 524,9102 652,7065 429,8729 710,6260 729,8898 779,5086 760,5005 [8] 746,7065 675,1872 809,3881 826,6930 702,2512 725,9460

[[11]] NULL

[[12]] [1] 92,74722 90,60442 84,34888 80,32372 142,12306 95,99942 78,94941 [8] 115,79662

[[13]] [1] 176,4084 382,2073 293,2750 217,3138 196,9911

[[14]] [1] 964,98822 497,56903 457,62909 807,84265 762,32831 171,38272 386,63550 [8] 97,36595 48,06548 26,26266

[[15]] [1] 95,6325 126,6646 141,5341 131,3787 105,1319 111,2801 188,6825 [8] 269,6385 188,2933 125,9783 133,8734 272,3000

[[16]] [1] 714,4169 371,5448 444,0142 626,7087 280,6295 547,9183 355,3295

[[17]]

228,8269

[[18]] [1] 114,6883 191,4677 153,4848 103,0313 412,7148 267,1815

[[19]] [1] 81,57622 75,48589

[[20]] [1] 831,36621 800,62671 16,76544 168,52444 397,27487 534,42041 62,76883 [8] 200,08041 56,51281

R specific loops

Ask for maximal value of elevation by county

lapply(dem_county, max)[1:4]

Warning in FUN(X[[i]], ...): no non-missing arguments to max; returning -
Inf

[[1]] [1] 773,7128

[[2]] [1] 356,1018

[[3]] [1] 144,0656

[[4]] [1] 387,1294

Example: static map - locally stored data

First we need to prepare data

#import geo data into R
library(rgdal)
#read via GDAL/OGR
#http://thematicmapping.org/downloads/world_borders.php
W <- readOGR(getwd(),"TM_WORLD_BORDERS-0.3", verbose = F, stringsAsFactors=FALSE) 
#or
#W <- download.file("http://www.mappinghacks.com/data/TM_WORLD_BORDERS-0.3",getwd(), method="auto",mode="wb")

Inspect the data

#inspect CRS/SRID
W@proj4string

CRS arguments: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

#take a look at the data
head(W@data[1:5], n=3)

FIPS ISO2 ISO3 UN NAME 0 AC AG ATG 28 Antigua and Barbuda 1 AG DZ DZA 12 Algeria 2 AJ AZ AZE 31 Azerbaijan

Base graphic system, static map

plot(W, main="Hello World") 

Base graphic system, improvement of static map

plot(W, main="Hello World", col='olivedrab3',bg='lightblue') 

Base graphics

Zoom, additional improvements

#create gridlines
grat <- gridlines(W, easts=seq(-40,20,by=5), norths=seq(10,80,by=5), ndiscr = 20)
#zoom to specified spatial extent on plot
plot(W, xlim=c(-30,30), ylim=c(30,58),col='olivedrab3',bg='lightblue') 
lines(grat, lty=2)#add gridlines
title(main="My map with zoom", sub="Gridlines added") #add title

Additional improvements

Base graphic, different adjustments

load("df.RData") #load data.frame saved in R format
#create spatial objects from data.frame
coordinates(df) <- ~ lng  + lat
#define CRS to the points
#World Geodetic System - WGS
proj4string(df) <- CRS("+proj=longlat +datum=WGS84") #EPSG:4326
#colors(2)
col_2 <- colors()[351]
x <-Polygon(cbind(x=c(13.9, 20.2, 20.2, 13.9), y=c(42.5, 42.5,46.5, 46.5)))
class(x)

[1] “Polygon” attr(,“package”) [1] “sp”

str(x)

Formal class ‘Polygon’ [package “sp”] with 5 slots ..@ labpt : num [1:2] 17 44,5 ..@ area : num 25,2 ..@ hole : logi TRUE ..@ ringDir: int -1 ..@ coords : num [1:5, 1:2] 13,9 20,2 20,2 13,9 13,9 42,5 42,5 46,5 46,5 42,5

Different map, base graphics

• basline data the same

plot(W,  col=col_2, xlim=c(13.9, 19.6), ylim=c(42.5, 46.5), main="Hello Croatia", lwd=2,  axes=T,  xlab="Longitude", ylab="Latitude", col='lightgray',bg='lightblue')
#add points to existing graph
points(df,  pch=21, cex=1.75, col="red4")
axis(1, col.axis = "dark gray")
axis(2, col.axis = "dark gray")
box(lty = 'solid', col = 'light gray')
text(x=14.7, y=46, labels="AUSTRIA", cex=0.75)       
text(x=17.8, y=44.5, labels="BOSNIA AND HERZEGOVINA", cex=0.75)   
text(x=18, y=46.4, labels="HUNGARY", cex=0.75) 
text(x=19.2, y=42.7, labels="MONTENEGRO", cex=0.75)

Improved map

Leaflet example, selected polygons/counties

• Oroslavje and Lobor

m <- leaflet()
m <- addTiles(m) %>%
addPolygons(data=W,color = "#03F", weight = 5, opacity = 0.5,) %>%
addPolygons(data=select,color = "#03F", weight = 1, opacity = 0.5,)
m

Example: data from MS Access

• Not spatial DB ..

#connect to the DB
ptice.db <-odbcDriverConnect("ptice_gps.mdb")
#SQL query
galeb_52 <- sqlQuery(ptice.db, query="SELECT  (latitude) AS lat, (longitude) AS lon , gps_fixstatus ,altitude, date_time    FROM galebovi_2008  WHERE serijski_broj = 52")
galeb_52_3d <- galeb_52[which(galeb_52$gps_fixstatus == '3D'),] 

Create ‘spacetime’ class object - package ‘spacetime’

sp <- SpatialPoints(galeb_52_3d[,c("lon","lat")])
proj4string(sp) <- CRS("+proj=longlat +datum=WGS84")
library(spacetime)
galeb_52_3d.st <- STIDF(sp, time = galeb_52_3d$date_time, data = galeb_52_3d)
#select and sort by time
galeb_52_3d.st_600   <- galeb_52_3d.st[1:600,]

#shepe of an icon
shape <- "R.png"

#export as kml
library(plotKML)

plotKML version 0.5-8 (2017-05-12)

URL: http://plotkml.r-forge.r-project.org/

kml(galeb_52_3d.st_600, labels = date_time, colour = altitude,shape=shape, colour_scale=rep("#FF00FF", 4), kmz = F)

KML file opened for writing...

Writing to KML...

Closing  galeb_52_3d.st_600.kml

Example: data from PostGres/PostGIS

library(RPostgreSQL) #R / PostgreSQL connect
library(RODBC) #connect to diff ODBC DB
library(rpostgis) #additional functionality for PostgreSQl/PostGIS
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, user="postgres",  dbname = "GEO_meetup", host = "PostgreSQL9.6(localhost:5432)", password="*****")
dbListTables(con)
HRV_adm1 <- "SELECT * FROM HRV_adm1*)
dbDisconnect(con)

Examples of analysis, points neighbours (sp)

Nearest neighbourgs, k=2; library(spdep)

library(spdep)
HRV_adm2_nn2<-knearneigh(coordinates(HRV_adm2), k=2) 

Visualisation of neighbours

plot(HRV_adm2, border="darkgrey", axes=T)
plot(knn2nb(HRV_adm2_nn2), coordinates(HRV_adm2), add=TRUE, col="red4")
title(main="K nearest, k = 2")

Analysis of point pattern

p_HR_300 <-  spsample(HRV_adm2,300, "random") 
class(p_HR_300)

[1] “SpatialPoints” attr(,“package”) [1] “sp”

points_ppp <- as.ppp(p_HR_300) #create point pattern
points_ppp

Planar point pattern: 300 points window: rectangle = [13,556343, 19,210248] x [42,53421, 46,50464] units

den_10000 <-  density(points_ppp, 10000) #density of point pattern at 10000 m
summary(den_10000)#summary of the result

real-valued pixel image 128 x 128 pixel array (ny, nx) enclosing rectangle: [13,55634, 19,21025] x [42,53421, 46,50464] units dimensions of each pixel: 0,0442 x 0,031019 units Image is defined on the full rectangular grid Frame area = 22,4484489889368 square units Pixel values range = [13,36395, 13,36395] integral = 300 mean = 13,36395

Visualisation of result, export is necessary

plot(density(points_ppp), main="Intensity of point process", col=gray(0:40/40))
points(p_HR_300,  cex=1.2, pch=19, col="darkred",add==T)

Library geosphere

When data are not in Cartesian system

pts <- SpatialPoints(data.frame(x = c(8.33, 30), y = c(47.22, -23)))
proj4string(pts) <- CRS("+proj=longlat +datum=WGS84")

# compute the great-circle
gcLine <- greatCircle(pts[1], pts[2], sp = TRUE)

Visualising results

#png("great-circle.png", width=10, height=10, units="cm", res=300 )
spplot(W[1], col.regions = "grey80", sp.layout = list(list("sp.lines", gcLine, 
    lwd = 2, col = "blue"), list("sp.points", pts, pch = 16, cex = 1.5, col = "red")), colorkey = FALSE, main="Great-circle trajectory")
#invisible(dev.off())
rm(img)
img <- readPNG("great-circle.png")
grid.raster(img)

Call external software for some algorithm

• Example calling SAGA GIS geoprocessor
  • Possible to have several versions of SAGA, specify in the call
• QGIS
• ESRI …

library(RSAGA)
rsaga.get.modules(env = rsaga.env() )
myenv <-  rsaga.env(path="C:\\Program Files (x86)\\SAGA-GIS")
rsaga.get.modules(env=myenv)
rsaga.get.modules("shapes_polygons")
rsaga.get.usage("shapes_polygons", 2)
rsaga.geoprocessor (lib="shapes_polygons", 2, param=list(POLYGONS="my_shapefile.shp", OUTPUT="result_shapefile.shp", BAREA=T))